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Abstract 

Poisson Voronoi diagrams are useful for modeling and describing various natural 
patterns and for generating random lattices. Although this particular space tes- 
sellation is intensively studied by mathematicians, in two- and three dimensional 
spaces there is no exact result known for the size-distribution of Voronoi cells. Mo- 
tivated by the simple form of the distribution function in the one-dimensional case, 
a simple and compact analytical formula is proposed for approximating the Voronoi 
cell's size distribution function in the practically important two- and three dimen- 
sional cases as well. Denoting the dimensionality of the space by d (d = 1, 2, 3) the 
f(y) = Const * y^ 3d ^ 1 '' 2 exp(—(3d + V)y/2) compact form is suggested for the nor- 
malized cell-size distribution function. By using large-scale computer simulations 
the validity of the proposed distribution function is studied and critically discussed. 

Key words: Voronoi diagrams, Monte Carlo methods, cell-size distribution 
PACS: 89.75.Kd, 02.50.Ng, 02.10.-v 



O 



o 
o 



S^ 1 Introduction 
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Voronoi diagrams [1] are a particular case of space tessellation, where given a 
set of centers, the space is divided according to their "spheres of influence". 
Each Voronoi cell contains those points of the space that are closest to the same 
center. A Voronoi tessellation in two-dimension would look like the polygons 
sketched in Figure 1 or Figure 2d. 
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Fig. 1. The "perpendicular bisector method" for constructing Voronoi diagrams in 
2D. 

Given a set of centers there are two relatively easy ways to generate the cor- 
responding Voronoi diagram. We sketch this methods for the two-dimensional 
(2D) case, and the generalization to any other dimension is immediate. In 
the first method (the perpendicular bisectors method [1,2]) one starts from a 
given center (Pq) and detects the nearest Pi center to it. A part of the per- 
pendicular bisector on the P$P\ line will form the first edge of the Voronoi 
polygon corresponding to Pq. Than the second nearest center (P2) is detected 
and the perpendicular bisector on P0P2 is constructed again. This algorithm 
is continued with the third (P3), fourth (P4), fifth (P5),... nearest center, until 
the perpendicular bisectors on P0P3, P0P4, P0P5.... will close a stable polygon 
which does not changes after considering any more distant points. Repeating 
the above algorithm for all centers the Voronoi tessellation of the whole space 
(Figure 1) can be obtained. 

The second method (called the Avrami- Johnson- Mehl method [3]) is especially 
useful for computer simulations. In this algorithm each center is identified as 
a nucleation point, from where a virtual disc with uniform radial velocity is 
growing (Figure 2). When two discs touch each other the growth in the contact 
direction is stopped for both of them and the contact point becomes a point 
on the corresponding Voronoi diagram. The growth in all other directions is 
continued until a nearby disc is reached. In this way the same space tessellation 
as in the perpendicular bisector algorithm is achieved. In computer simulations 
it is handy to implement this algorithm not on continuous space but on large 
lattices since the contact points are easier to identify. 



In the two-dimensional case the Voronoi diagram can be obtained also from 




Fig. 2. The Avrami-Johnson-Mehl method for constructing Voronoi diagrams in 2D. 
Figures a.-d. presents snapshots from a small graphical simulation. 

Delaunay triangulation. The Delaunay triangulation of a point set is a collec- 
tion of edges satisfying an empty-circle criteria, which means that for each edge 
we can find a circle containing the edge's endpoints but not containing any 
other point from the initial set. In two-dimension the Delaunay triangulation 
is the dual structure of the Voronoi diagram [2]. 

A particular case of Voronoi diagrams, where the centers are randomly and 
uncorrelated distributed, are called Poisson Voronoi diagrams. Poisson Voronoi 
Diagrams (PVD) are especially important for modeling and describing a wide 
variety of natural and social phenomena. PVD has been used to construct 
random lattices in quantum field theory [4] or in the studies of conductivity 
and percolation in granular composites [5]. PVD was also used in modeling 
growth of metal clusters on amorphous substrates [6], in studying conduction 
and percolation in continuous media [7], in modeling microemulsions [8], in 
interpreting small angle X-ray scattering for heterogeneous catalyst [9], in 
evaluating the actual galaxy distribution [10], in describing sections through 
various geological materials [11], in biology [12], in animal ecology [13], in 
sociology [14] etc... The above list is far from being complete, and suggests 
just a few possible applications for this particular space tessellation. For a 
more complete discussion of the use of Voronoi diagrams many good review 
works are available [15,16,17]. 

Despite their importance in science, our knowledge on the geometrical and 
statistical properties of PVD is far from being complete [1,2]. One of the 
most debated and less clarified aspect is the g(S) size distribution function of 
Voronoi cells g(S) = P(S, S + dS)/dS, where P(S, S + dS) is the probability 



that the size of a Voronoi cell is between S and S+dS. Instead of g(S) it is more 
convenient to use the more general f(y) distribution function for the y = S/(S) 
normalized cell sizes, which is independent of the center's density and it is 
universal for all PVD in a given dimension. Alternatively, one could determine 
the F(y) cumulative distribution function defined as F(y) = Jq f(x)dx. 

Apart of the simple one-dimensional case, presently there is no exact result or 
handy analytical approximation for the form of f(y). Since f(y) in two and 
three dimensions is of primary interest in many practical applications, there is 
a growing need for a simple and analytically usable formula. This would help 
characterizing and classifying several experimental patterns, and would give an 
important starting point also for modeling these structures. In contrast with 
mathematicians experimental scientist need a simple expression that could 
give a first hint about the nature of the measured cell-size distribution, which 
is usually determined with a poor statistics. 

There are many conjectures on the analytical form of f(y) and many computer 
simulations were done to prove the suggested forms. Up to our knowledge in 
2D the largest computer simulations were done by Tanemura [18,19] with 10 7 
Voronoi cells and Hinde and Miles [20] with 2 x 10 6 cells. In 3D the largest 
ensembles were studied again by Tanemura [18,19] (3 x 10 6 cells) and Kumar 
et al. [21] (3.6 x 10 6 cells). 

As a generally accepted result emerges, that a three parameter (a, b and c) 
generalized gamma function fit 

f(y) = c-^y^expi-byc), (1) 

T(a/c) 

describes the computer simulation data reasonable well. Some authors [22,6] 
suggested however that a simpler two-parameter (a and b) gamma function fit 

f(y) = SrV^expi-by) (2) 

r(o) 

works also well. 

In 2D for the three parameter fit (1) Tanemura [18,19] found a = 3.315, 
b = 3.04011 and c = 1.078, in good agreement with the results of Hinde and 
Miles [20] a = 3.3095, b = 3.0328 and c = 1.0787. For the two parameter fit 
(2) the values a = b = 3.61 [22] or a = 3.61 and b = 3.57 [6] were reported. 

In the 3D case Tanemura found [19] a = 4.8065, b = 4.06342 and c = 1.16391 
for (1), while Kiang [23] suggested a fit of the form (2) with a = b = 6. We 
have to mention however that the simulations of Tanemura [19,18] did not 
supported Kiang's results [23] at all. 



For the sake of completeness it also has to be mentioned that there is an 
exact analytical result for the second moment ({y 2 )) of the PVD's both in 
the two and three-dimensional cases [24]. According to this (y^o) = 1-280 and 
(y 2 D ) = 1.180 [15], offering an excellent possibility for testing the computer 
simulation results and the correctness of the proposed fit. It was the enormous 
discrepancy between Gilbert's and Kiang's results for this second moment that 
condemned Kiang's simulation results. 

The aim of the present work is not to give a better and more complicated fit 
for f(y). We would rather intend to prove that a simple two parameter fit of 
the form (2) used by Kiang can be still a fair approximation for all practical 
applications. Experimental scientist instead of focusing on a more accurate 
but difficult fit for the presumed PVD type patterns can use with confidence a 
simple approximation of the form (2). In the present work large-scale computer 
simulations are also considered for the problem, generating more Voronoi cells 
than in all previous studies we are aware of. The statistics of 3 x 10 7 and 
1.8 x 10 7 cells are studied in 2D and 3D, respectively. Using this improved 
statistics Kiang's conjecture will be followed and a first approximation for 
f(y) in the form (2) will be given with simple and handy values of a = b. 
Solving the problem exactly in ID will give us further motivation for this 
simpler form of f(y). 



2 The one-dimensional case 



Let us first study theoretically the simple problem in ID and prove the validity 
of (2) with a = b. One has to mention however that several other methods 
are known to obtain the exact form of the fmiy) distribution function in this 
simple case [3,15]. 

A line with length L is considered on which randomly and independently N 
centers are distributed. The density of centers is given thus as n = ^- = -X-, 
where (d) stands for the average distance between centers. We will study the 
limit L — > oo, N — * oo, but n finite. It is immediate to construct the Voronoi 
diagrams for these centers (Figure 3). If a center P is considered, first it's 
neighbor in the left (Pi) and right (P r ) direction will be detected. Than the 
PPi and PP r lines are divided in two equal parts, by the D\ and D r points, 
respectively. The segment D[D r is than the Voronoi cell corresponding to the 
center P. It is obvious that for the considered limit the average length of 
Voronoi cells is (d) = 1/n. 

In order to get the distribution function g(d) of the Voronoi cell's length, 
first the distribution function h(s) for the lengths between centers will be 
determined. Let us start from the well-known Poisson distribution P(N,t), 
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Fig. 3. Construction of Voronoi cells in ID. 
giving the probability that inside a length t there are N centers: 

P(N,t) = ±{N)?exp((N) t ). (3) 

In the above equation (N) t = nt stands for the expected (average) number of 
centers on a length t. The probability that in an interval of length t situated 
on the immediate right of P there are no other centers is: 

P(0,t) = exp(-nt). (4) 

The cumulative distribution P r (d r > t) that the first neighbor at the right is 
at a distance d r bigger than t is P r (d r > t) = P(0, t). The distribution function 
g r (d r ) for the lengths d r can be thus calculated as: 

g r (d) = r —^- = n exp(-nd). (5) 

Due to symmetry constraints, the same distribution function should apply for 
the di lengths relative to the first neighbor in the left direction. The distri- 
bution function for the half of these intervals (z = d r /2 or z = di/2) is given 

as: 

w(z) = 2nexp(—2nz). (6) 

The length d of the Voronoi cell is d = -j + y , and it's distribution g(d) can 
be calculated as the convolution of two distributions of form (6): 

g{d) = I w(z)w{d — z)dz = Anexp{—2nd). (7) 

Jo 



It is immediate to realize that this distribution function is normalized for 
L — ► oo. The distribution function for the adimensional quantity y = d/(d) is 
given then as 

f 1D (y) = 4yexp(-2y), (8) 

which has the general form (2) with a = b = 2. The cumulative distribution 
function FiE>(y) is given by 

F 1D (y) = l-{2y + l)e- 2u , (9) 



and the moments of fmiy) are immediately calculable: {y)m = 1; (y 2 )iD 

3. 
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:: " (y 3 )iD — 3. The most probable normalized length obtained from (8) is 




Fig. 4. Simulation results (empty circles) in ID in comparison with the (8) exact 
result (solid line). Results for the cumulative distribution function are also plotted. 
Filled circles are simulation data and the dashed line is given by equation (9). 

A simple computer simulation exercise can easily convince us about the valid- 
ity of our calculations. Results in this sense are presented on Figure 4. As an 
interesting observation one can realize that the distribution function for the 
lengths between randomly displaced centers (given by (5) ) is qualitatively 
different from the (7) distribution function for the length of Voronoi cells (see 
also [25]). 



3 The two-dimensional case 



Theoretical attempts to get analytical result for f 2 D (y) (y = S/(S), with S 
the area of Voronoi cells) in 2D, failed. We considered thus Monte Carlo-type 
computer simulations and fitted our simulation data in different forms. In 
particular, we focused on a three parameter fit in the generally accepted (1) 
form and tried also a simple two parameter approximation (2) with handy 
a = b values. It was found that the simple choice a = b = 7/2 gives a visually 
good fit. For the normalized distribution function of Voronoi cell areas in 2D 
we proposed thus the 



hn{y) = ^§^V**°P(-lv) 



(10) 



simple approximation. On Figure 5 we plotted with a continuous line the curve 
(10) in comparison with simulation data obtained on 29.889 x 10 6 Voronoi cells 
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Fig. 5. The Voronoi cell's normalized area-distribution function in 2D. Empty circles 
are simulation results and the solid line is the (10) formula suggested in this study. 
With point-dashed line the cumulative distribution function, and with dashed line 
the (1) best gamma function fit is plotted. The two inset figures are magnification 
of the small and large y limits plotted on log-log and log-normal scales, respectively. 
On the scale of the figure there is no detectable difference between the cumulative 
distribution function calculated from simulation and the analytical expression given 
by (1) and (10). 

(almost three times more than the number of cells used by Tanemura). On 
the same graph it is drawn with dashed line the best gamma function (1) fit. 
The F2d{d) cumulative distribution function is also plotted with dash-dotted 
line. 



At a first glance there is no detectable difference between computer simulation 
results, the curve suggested by (10) and the gamma-function fit. Magnifying 
however the initial part and tail of the distribution function and plotting 
it on log-log and log-normal scales (insets in fig. 5), respectively, one can 
observe slight differences. As expected, the three parameter gamma-function 
fit is better, but the improvement relative to (10) is not spectacular. The 
best-fit parameters obtained by us for (1) are a = 2.2975, b = 3.01116 and c = 
1.0825, in comparison with the values a = 3.315, b = 3.04011 and c = 1.078 
obtained by Tanemura [19,18]. For the analytically known second moment 
of the distribution ((y^ eor2 z)) = 1-280) our simulation data gives (y 2 sim 2D) = 
1.28231, and the three-parameter gamma-fit yields {y^ amma 2D) = 1-27947. The 



error relative to the exactly known result is of the same order (0.04%) as in 
the case of the fit given by Tanemura to his own computer simulation results. 

Using (10) all the important moments can be analytically calculated: (y)2D = 
1; (y 2 )2D = f ; (y 3 )2D = f§- The second moment has of course a much bigger 
relative error (0.4%) respective to the exactly known value than the one ob- 
tained with the more sophisticated three-parameter gamma-function fit. This 
relative error is however still quite small and usual experimental data on Pois- 
son Voronoi type patterns gives deviations of the order of a few percents. The 
most probable normalized area is y pro b2D = f • 



4 The three-dimensional case 



Due to the complex geometry involved, the possibility to analytically calcu- 
late fsDiy) {y = V/(V), V the volume of Voronoi cells) in 3D is even more 
gloomy. We performed thus again large scale computer simulations, studying 
the statistics of 18.27 x 10 6 Voronoi cells (six times more than the statistics 
considered by Tanemura). The three-parameter gamma- function gives a good 
fit for the simulation data, but again as in the two-dimensional case a simple 
fit of the form (2) works also reasonably well and handy a = b = 5 values 
can be considered. We suggest thus that in 3D the Voronoi cell's normalized 
volume distribution can be approximated as: 



fMv) = ~^y 4 exp(-5y) (11) 

On Figure 6 simulation results (empty circles) are compared with the (11) 
approximation (continuous line) and the three parameter gamma-function fit 
(dashed-line) . In a first visual approximation one will realize that both curves 
describe well the simulation data. Magnifying however the initial part and 
tail of the distribution function and plotting it on log-log and log-normal 
scales (insets in Figure 6), respectively, one can observe the differences. As 
expected, the (1) gamma- function fit is better, and follows more the trend 
of the simulation results. The best fit parameters obtained in this study are 
a = 3.24174, b = 3.24269 and c = 1.26861 (in contrast with a = 4.8065, 
b = 4.06342 and c = 1.16391 found by Tanemura [19,18]). The improvement 
relative to the simple (11) approximation is however again not spectacular, 
and is relevant only in the limit of very large or very small Voronoi cells. These 
limits does not appear usually in real experimental data, due to the fact that 
a much weaker statistics is achieved (patterns with less than 10 4 cells are 
studied). On the figure it is also plotted (point-dashed line) the form of the 
cumulative distribution function F^^y). On the scale of the image there is no 
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Fig. 6. The Voronoi cell's normalized volume-distribution function in 3D. Empty 
circles are simulation results and the solid line is the (11) formula suggested in 
this study. With point-dashed line the cumulative distribution function, and with 
dashed line the (1) best gamma function fit is plotted. The two inset figures are 
magnification of the small and large y limits plotted on log-log and log-normal scales, 
respectively. On the scale of the figure there is no detectable difference between 
the cumulative distribution function calculated from simulation and the analytical 
expression given by (1) and (11). 



detectable difference in the cumulative distribution function determined from 
simulation and the forms (11) or (1). 



By using 

i; (y 2 ) 



3D 



'11) the important moments are analytically calculable: (y)$D 
{y 3 )3D = §§• The most probable normalized volume is y pro b3D 
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|. For the second moment the relative error respective to the analytical exact 



results ((ytheor3D = 1-18) is e = 1.7%. The gamma-fit for the simulation data 
yields {ylogamma) = 1-18683 (e = 0.57%) while Tanemura's fit seems better 
yielding 1.17830 (e = 0.14%). The second moment computed directly from 
simulation data is (y^ im3 D) = 119, giving the e = 0.85% relative error. 



In agreement with the simulations of Tanemura [19,18] we have also found 
that the values a = b = 6 suggested by Kiang [23] are not appropriate and 
give no good fit to our simulation data. 
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5 Conclusions 



Motivated by the simple form of the exact result (8) for the size distribu- 
tion function of Poisson Voronoi cells in ID we proposed simple expressions 
for approximating the distribution in 2D and 3D where no exact results are 
available. Exceeding the statistics considered in all previous studies computer 
simulations were used to investigate numerically the distribution function. It 
was shown that a simple form (2) with a = b is appropriate for all practical 
applications to approximate the size-distribution of the Poisson Voronoi cells. 
In ID the exact results gives a = 2. In 2D and 3D we found that a = 7/2 
and a = 5, respectively, gives fair approximation. The simple values suggested 
for a = b allows also to write in a compact form the approximations (10) and 
(11). If we denote by d the dimensionality of the problem (d = 1,2,3), the 
value of a can be given as a = (3d + l)/2. Equations (8), (10) and (11) can be 
written then in a compact form as: 

„ /x ((3rf+ i)/ 2 )(3^+i)/2 3d _ x 3d+1 

fd(y) = - — „ , \ , ^ — y—exp( y) (12) 

Jyy> r((3d+l)/2) y v 2 yj K J 

This distribution function is not an exact one and it is less accurate than 
a more complicated three parameter fit given by the generalized gamma- 
function. Mathematicians will probably not appreciate it... but due to it's 
simplicity it will definitely be of importance for experimental scientist study- 
ing and characterizing complex Voronoi diagram-like patterns. 
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